%% initialization
ROOTDIR = '..';
MEXBINDIR = [ROOTDIR filesep 'bin'];
[~,~,~] = mkdir(MEXBINDIR);
addpath([ROOTDIR filesep 'common'], MEXBINDIR, [ROOTDIR filesep 'IBP'], fullfile(ROOTDIR, '..', 'daSVM'));

N = 10;
alphav = 4;

maxiter = 200;
numiter = 500;

gfuncs = {
    @(Z)size(Z,2), ...
    @(Z)mean(sum(Z,2)), ...
    @(Z)sum(sum(Z) == 1), ...
    @(Z)sum(sum(Z) == (N+mod(N,2))/2), ...
    @(Z)numel(unique(2.^(N-1:-1:0)*Z)) };

gfuncnames = {
    '# of features', ...
    'avg. # of active features by each object', ...
    '# of singleton features', ...
    '# of half-full features', ...
    '# of unique histories'};

gavgvals = [alphav*hrmsum(N), alphav, alphav, alphav*2/(N+mod(N,2)), nan];
gstdvals = [sqrt(alphav*hrmsum(N)), nan, sqrt(alphav), sqrt(alphav*2/(N+mod(N,2))), nan];

mrgcndstats = zeros(numiter, numel(gfuncs));
succndstats = zeros(size(mrgcndstats));
avgmrgcndstats = zeros(size(mrgcndstats));
stdmrgcndstats = zeros(size(mrgcndstats));
avgsuccndstats = zeros(size(succndstats));
stdsuccndstats = zeros(size(succndstats));

savedir = input('save to: ', 's');
mkdir(savedir);

%% getting-it-right
Z2 = ibprnd(N, alphav);
for i = 1:numiter
    Z1 = ibprnd(N, alphav);
    mrgcndstats(i,:) = cellfun(@(f)f(Z1), gfuncs);
    
    for j = 1:maxiter
        Z2 = ibpcrnd(Z2,alphav,randperm(N));
    end
%     Z2 = ibphrnd(N, alphav);
    succndstats(i,:) = cellfun(@(f)f(Z2), gfuncs);
    
    avgmrgcndstats(i,:) = mean(mrgcndstats(1:i,:),1);
    stdmrgcndstats(i,:) = std(mrgcndstats(1:i,:),1,1);
    avgsuccndstats(i,:) = mean(succndstats(1:i,:),1);
    stdsuccndstats(i,:) = std(succndstats(1:i,:),1,1);
    
    stats = [
        num2str([avgmrgcndstats(i,:)', stdmrgcndstats(i,:)'], ...
                ['%.3f',char(177),'%.3f']), ...
        repmat(blanks(4), numel(gfuncs), 1), ...
        num2str([avgsuccndstats(i,:)', stdsuccndstats(i,:)'], ...
                ['%.3f',char(177),'%.3f'])
    ];
    fprintf('%d:\n', i);
    disp(stats);
end

save([savedir filesep 'stats.mat'], 'N', 'alphav', ...
    'maxiter', 'gfuncs', 'gfuncnames', 'gavgvals', 'gstdvals', 'numiter', ...
    'mrgcndstats', 'avgmrgcndstats', 'stdmrgcndstats', ...
    'succndstats', 'avgsuccndstats', 'stdsuccndstats');

%% display the results
for j = 1:numel(gfuncs)
    if isnan(gavgvals(j))
        figure;
        plot(avgmrgcndstats(:,j)+stdmrgcndstats(:,j), '-.r'); hold on;
        plot(avgmrgcndstats(:,j)-stdmrgcndstats(:,j), '-.r');
        plot(avgmrgcndstats(:,j), '-.b');
        plot(avgsuccndstats(:,j)+stdsuccndstats(:,j),'r');
        plot(avgsuccndstats(:,j)-stdsuccndstats(:,j),'r');
        plot(avgsuccndstats(:,j),'b'); hold off;
        title([gfuncnames{j}, ': marginal/successive -conditional']);
        pause
        continue;
    end
    subplot(2,1,1);
    plot(avgmrgcndstats(:,j)+stdmrgcndstats(:,j), 'r'); hold on;
    plot(avgmrgcndstats(:,j)-stdmrgcndstats(:,j), 'r');
    plot(avgmrgcndstats(:,j), 'b');
    plot(1:size(avgmrgcndstats,1), gavgvals(j)+gstdvals(j), '-.r');
    plot(1:size(avgmrgcndstats,1), gavgvals(j)-gstdvals(j), '-.r');
    plot(1:size(avgmrgcndstats,1), gavgvals(j), '-.b'); hold off;
    title([gfuncnames{j}, ': marginal-conditional']);
    
    subplot(2,1,2);
    plot(avgsuccndstats(:,j)+stdsuccndstats(:,j),'r'); hold on;
    plot(avgsuccndstats(:,j)-stdsuccndstats(:,j),'r');
    plot(avgsuccndstats(:,j),'b');
    plot(1:size(avgsuccndstats,1), gavgvals(j)+gstdvals(j), '-.r');
    plot(1:size(avgsuccndstats,1), gavgvals(j)-gstdvals(j), '-.r');
    plot(1:size(avgsuccndstats,1), gavgvals(j), '-.b'); hold off;
    title([gfuncnames{j}, ': successive-conditional']);
    pause
end
close all